// @HEADER
//*********************************************************************//
//  SiYuan: A numerical PDE solver                                     //
//  Copyright (2022) YUAN Xi                                           //
//  This Software is released under the BSD 2-Clause license detailed  //
//  in the file "LICENSE" in the top-level SiYuan directory            //
//*********************************************************************//
// @HEADER

#ifndef _GENEIGENSOLVER_H
#define _GENEIGENSOLVER_H

#include <iostream>

//#include "LOCA.H"
#include "Thyra_ResponseOnlyModelEvaluatorBase.hpp"
#include "Thyra_ModelEvaluatorDefaultBase.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_ParameterList.hpp"

#include "Tpetra_Operator.hpp"
#include "Amesos2.hpp"
#include "AnasaziOperator.hpp"

#include "EigenObserver_WriteToExodus.hpp"

namespace Anasazi {

template<class Scalar>
class TpetraGenOp : public virtual Tpetra::Operator<Scalar> {

  typedef Tpetra::CrsMatrix<Scalar>      CrsMatrix;
  typedef Tpetra::MultiVector<Scalar>    MV;

public:
  // \brief Constructor
  //
  // \param solver [in] The Amesos2 solver to use for solving linear
  //   systems with K.
  // \param massMtx [in] The "mass matrix" M.
  // \param useTranspose [in] If false, apply $K^{-1} M$; otherwise,
  //   apply $M K^{-T}$.
  TpetraGenOp (const Teuchos::RCP<Amesos2::Solver<CrsMatrix,MV>>& solver,
               const Teuchos::RCP<CrsMatrix>& massMtx,
               const bool useTranspose = false );

  // The Operator's (human-readable) label.
  const char* describe() const {
    return "Operator that applies K^{-1} M or M K^{-T}";
  }

  Teuchos::RCP< const Tpetra::Map<> > getDomainMap() const {
    return massMtx_->getDomainMap();
  }

  Teuchos::RCP< const Tpetra::Map<> > getRangeMap() const {
    return massMtx_->getRangeMap();
  }

  //! Apply method [inherited from Tpetra::Operator class]
  // \brief Compute Y such that $K^{-1} M X = Y$ or $M K^{-T} X = Y$
  // , where solves with the Tpetra::CrsMatrix K use an Amesos
  //   solver.
  void apply(const Tpetra::MultiVector<Scalar> &X, Tpetra::MultiVector<Scalar> &Y,
           Teuchos::ETransp mode = Teuchos::NO_TRANS,
           Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
           Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const override;

private:
  // Default constructor: You may not call this.
  TpetraGenOp() {};

  Teuchos::RCP< Amesos2::Solver<CrsMatrix,MV> > solver_;
  Teuchos::RCP<CrsMatrix> massMtx_;
  bool useTranspose_;
};

//! Partial specialization of OperatorTraits for TpetraGenOp objects.
template <class Scalar, class LO, class GO, class Node>
class OperatorTraits<Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node>, TpetraGenOp<Scalar> >
{
  public:
    static void
    Apply (const TpetraGenOp<Scalar>& Op,
           const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
           Tpetra::MultiVector<Scalar,LO,GO,Node>& Y)
    {
      Op.apply(X, Y, Teuchos::NO_TRANS);
    }
};

}  // namespace Anasazi



namespace SiYuan {

template<class Scalar> class EigenSolver;

/** 
 *   Pivot value
*/
struct S_Pivot
{
  double pk;

  S_Pivot(double p) : pk(p) {};
};


/** \brief Nonmember constuctor.
 *
 * \relates EigenSolver
 */
template<class Scalar>
Teuchos::RCP<EigenSolver<Scalar> >
buildEigenSolver( const Teuchos::RCP<Teuchos::ParameterList>& pl,
  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > &model ,
  const std::shared_ptr< SiYuan::EigenObserver_WriteToExodus<Scalar> >&);

/** \brief Generalized Eigen ModelEvaluator.
 *
 * This class is used to calcualte Eigen values
 */
template<class Scalar>
class EigenSolver : public Thyra::ResponseOnlyModelEvaluatorBase<Scalar>
{
  typedef Thyra::ModelEvaluatorBase MEB;

public:

  /** \name Initializers/Accessors */
  //@{

  /** \brief . */
  EigenSolver( const Teuchos::RCP<Teuchos::ParameterList>& pl,
    const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > &model,
    const std::shared_ptr< EigenObserver_WriteToExodus<Scalar> >&);

  //@}

  /** \name Public functions overridden from ModelEvaulator. */
  //@{

  /** \brief . */
  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_space(int l) const;
  /** \brief . */
  Teuchos::RCP<const Teuchos::Array<std::string> > get_p_names(int l) const;
  /** \brief . */
  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_g_space(int j) const;
  /** \brief . */
  Teuchos::ArrayView<const std::string> get_g_names(int j) const;
  /** \brief . */
  Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
  //@}


private: // functions

  /** \name Private functions overridden from ModelEvaulatorDefaultBase. */
  //@{

  /** \brief . */
  Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
  /** \brief . */
  void evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
    const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
    ) const;

  //@}

private: // data members
  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > x_space_;
  
  Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > model_;
  std::shared_ptr< EigenObserver_WriteToExodus<Scalar> > observer_;
  int num_p_, num_g_;

  //Eigensolver parameters
  std::string method;
  bool bHermitian;
  std::string which;
  int nev, blockSize, maxIters;
  double conv_tol;
  bool verbose;
  bool debug;
  double sigma;  //shift-invert
  int nrestart;

  double K_pivot;
};


} // namespace SiYuan

#include "Eigensolver_impl.hpp"

#endif
